#### Claassen, C. & McLaren, L. (2021). Does Immigration Produce a Public Backlash or Public Acceptance? 
#### Time-Series, Cross-Sectional Evidence from 30 European Democracies. *British Journal of Political Science*
#### Code to replicate the analyses in the article 

library(texreg)
library(dplyr)
library(tidyr)
library(ggplot2)
library(RColorBrewer)
library(psych)
library(plm)
library(MASS)

# working directory
WD = rstudioapi::getActiveDocumentContext()$path 
setwd(dirname(WD))
print( getwd() )

# load and check data
migrat = read.csv("immigration opinion tscs dataset replication.csv")
names(migrat)
migrat$Country = as.factor(migrat$Country)
levels(migrat$Country)
table(migrat$Year)

# create panel dataset
mig.plm = pdata.frame(migrat, index = c("Country", "Year"))
sort(rowSums(table(index(mig.plm), useNA = "ifany")))
pdim(mig.plm) # n=30, T=39, N=1170

# create western europe dataset
cnts.drop = c("Croatia", "Czechia", "Hungary", "Poland", "Romania", "Slovakia", "Slovenia", "Latvia", "Estonia", "Lithuania")
migrat.weur = migrat[!migrat$Country %in% cnts.drop,]
migrat.weur$Country = as.factor(as.character(migrat.weur$Country))
mig.plm.we = pdata.frame(migrat.weur, index = c("Country", "Year"))
sort(rowSums(table(index(mig.plm.we), useNA = "ifany")))
pdim(mig.plm.we) # n=20, T=39, N=780

# create within and between measures
migrat$Imm_sup_btw = ave(migrat$Imm_sup_z, migrat$Country, FUN=function(x) mean(x, na.rm=TRUE))
migrat$Imm_sup_within = migrat$Imm_sup_z - migrat$Imm_sup_btw
migrat$Imm_sal_btw = ave(migrat$Imm_sal_z, migrat$Country, FUN=function(x) mean(x, na.rm=TRUE))
migrat$Imm_sal_within = migrat$Imm_sal_z - migrat$Imm_sup_btw
migrat$Noncit_stck_btw = ave(migrat$Noncit_stck_pc, migrat$Country, FUN=function(x) mean(x, na.rm=TRUE))
migrat$Noncit_stck_within = migrat$Noncit_stck_pc - migrat$Noncit_stck_btw
migrat$Immig_pc_btw = ave(migrat$Immig_pc, migrat$Country, FUN=function(x) mean(x, na.rm=TRUE))
migrat$Immig_pc_within = migrat$Immig_pc - migrat$Immig_pc_btw
migrat$Net_migr_btw = ave(migrat$Net_migr_pc, migrat$Country, FUN=function(x) mean(x, na.rm=TRUE))
migrat$Net_migr_within = migrat$Net_migr_pc - migrat$Net_migr_btw



#### Figure 1

col1 = rgb(1, 0.55, 0, 1)
col2 = rgb(0, 0.45, 0.7, 1)
cnts = levels(migrat$Country)

pdf("fig1.pdf", height=7.5, width=6.5)
par(mfrow=c(6,5), mar=c(1.2,1,0.2,0.2), tcl=-0.15, cex=0.9)
for(i in 1:length(cnts)) {
  plot(x=migrat[migrat$Country==cnts[i], "Year"],
       y=migrat[migrat$Country==cnts[i], "Imm_sup_z"],
       type="n", xlim=c(1988, 2017), ylim=c(-3, 3), xlab="", ylab="", main="", axes=FALSE)
  grid(col=rgb(0,0,0,0.5), lwd=0.5)
  text(cnts[i], x=1988, y=2.5, cex=0.75, adj=0)
  lines(x=migrat[migrat$Country==cnts[i], "Year"],
        y=migrat[migrat$Country==cnts[i], "Imm_sup_z"],
        col=col1, lwd=2)
  lines(x=migrat[migrat$Country==cnts[i], "Year"],
        y=migrat[migrat$Country==cnts[i], "Imm_sal_z"],
        col=col2, lwd=2)
  axis(side=1, at=c(1990,2000,2010), labels=c(1990,2000,2010), mgp=c(1,-0.2,0), cex.axis=0.65)
  axis(side=2, at=c(-2,0,2), mgp=c(1,0.1,0), cex.axis=0.65)
  box()
}
text(x=1989, y=migrat[migrat$Country==cnts[length(cnts)] & migrat$Year==1995, "Imm_sup_z"]-0.5, 
     labels="Immigration mood", adj=0, col=col1, cex=0.65)
text(x=1993, y=migrat[migrat$Country==cnts[length(cnts)] & migrat$Year==2005, "Imm_sal_z"]-0.9, 
     labels="Immigration concern", adj=0, col=col2, cex=0.65)
dev.off()



#### Models

## Table 1

form0 = diff(Imm_sup_z, lag=1) ~ plm::lag(Imm_sup_z, 1:2) + plm::lag(GDP_cp_grth, 1) + plm::lag(Unempl_WB, 1) 
form1 = update.formula(form0, . ~ . + ( diff(Net_migr_pc, lag=1) + plm::lag(Net_migr_pc, 1) ) * plm::lag(Noncit_stck_pc, 1))
form2 = update.formula(form1, . ~ . + plm::lag(FRight_seat_pr, 1) + plm::lag(BMI_i_ext, 1) + plm::lag(BMI_e_ext, 1) )
form3 = update.formula(form0, . ~ . + ( diff(Immig_pc, lag=1) + plm::lag(Immig_pc, 1) ) * plm::lag(Noncit_stck_pc, 1))
form4 = update.formula(form3, . ~ . + plm::lag(FRight_seat_pr, 1) + plm::lag(BMI_i_ext, 1) + plm::lag(BMI_e_ext, 1) )

perc.mods = lapply(list(form2, form4), plm, data=mig.plm, model="within", effect="indiv")
sapply(perc.mods, function(x) round(pwartest(x, order=1)[["p.value"]][[1]], 3) )
sapply(perc.mods, function(x) round(pcdtest(x, test="cd")[["p.value"]][[1]], 3) )
sapply(perc.mods, function(x) round(sigma(x), 3) )
vcm = lapply(perc.mods, function(x) plm::vcovSCC(x, cluster="time")) # cross-sectional heteroskedasticity and correlation
rse = lapply(vcm, function(x) sqrt(diag(x)))
pval = list()
for(i in 1:length(perc.mods)) {pval[[i]] = pt(abs(coef(perc.mods[[i]]) / rse[[i]]), df=500, lower.tail=FALSE) * 2}
screenreg(perc.mods, override.se=rse, override.pvalues=pval, digits=3, single.row=TRUE, leading.zero=FALSE, 
          stars=c(0.05))
vars = c("Immigration mood t-1", "Immigration mood t-2", "GDP growth per capita t-1", "Unemployment rate t-1",
         "Diff Net migration", "Net migration, % of pop. t-1", "Non-citizen stock, % of pop t-1",   
         "Far right seat share t-1", "Immigrant integration policy index t-1", "Immigrant entry policy index t-1",
         "Diff Net migration X non-citizen stock", "Net migration X non-citizen stock", 
         "Diff Immigration inflows", "Immigration inflows, % of pop. t-1", 
         "Diff Immigration inflows X non-citizen stock", "Immigration inflows X non-citizen stock")
htmlreg(perc.mods, file="table1.htm", override.se=rse, override.pvalues=pval, digits=3, 
        single.row=TRUE, leading.zero=FALSE, stars=0.05, custom.coef.names=vars,
        caption="Table 1: Drivers of change in immigration mood, 1-way FE mods, DKSEs")

## Table 2

form0 = diff(Imm_sal_z, lag=1) ~ plm::lag(Imm_sal_z, 1:2) + plm::lag(GDP_cp_grth, 1) + plm::lag(Unempl_WB, 1) 
form1 = update.formula(form0, . ~ . + ( diff(Net_migr_pc, lag=1) + plm::lag(Net_migr_pc, 1) ) * plm::lag(Noncit_stck_pc, 1))
form2 = update.formula(form1, . ~ . + plm::lag(FRight_seat_pr, 1) + plm::lag(BMI_i_ext, 1) + plm::lag(BMI_e_ext, 1) )
form3 = update.formula(form0, . ~ . + ( diff(Immig_pc, lag=1) + plm::lag(Immig_pc, 1) ) * plm::lag(Noncit_stck_pc, 1))
form4 = update.formula(form3, . ~ . + plm::lag(FRight_seat_pr, 1) + plm::lag(BMI_i_ext, 1) + plm::lag(BMI_e_ext, 1) )

sal.mods = lapply(list(form2, form4), plm, data=mig.plm, model="within", effect="indiv")
sapply(sal.mods, function(x) round(pwartest(x, order=1)[["p.value"]][[1]], 3) )
sapply(sal.mods, function(x) round(pcdtest(x, test="cd")[["p.value"]][[1]], 3) )
sapply(sal.mods, function(x) round(sigma(x), 3) )
vcm = lapply(sal.mods, function(x) plm::vcovSCC(x, cluster="time")) # cross-sectional heteroskedasticity and correlation
rse = lapply(vcm, function(x) sqrt(diag(x)))
pval = list()
for(i in 1:length(sal.mods)) {pval[[i]] = pt(abs(coef(sal.mods[[i]]) / rse[[i]]), df=500, lower.tail=FALSE) * 2}
screenreg(sal.mods, override.se=rse, override.pvalues=pval, digits=3, single.row=TRUE, leading.zero=FALSE, 
          stars=c(0.05))
vars = c("Immigration concern t-1", "Immigration concern t-2", "GDP growth per capita t-1", "Unemployment rate t-1",
         "Diff Net migration", "Net migration, % of pop. t-1", "Non-citizen stock, % of pop t-1",   
         "Far right seat share t-1", "Immigrant integration policy index t-1", "Immigrant entry policy index t-1",
         "Diff Net migration X non-citizen stock", "Net migration X non-citizen stock", 
         "Diff Immigration inflows", "Immigration inflows, % of pop. t-1", 
         "Diff Immigration inflows X non-citizen stock", "Immigration inflows X non-citizen stock")
htmlreg(sal.mods, file="table2.htm", override.se=rse, override.pvalues=pval, digits=3, 
        single.row=TRUE, leading.zero=FALSE, stars=0.05, custom.coef.names=vars,
        caption="Table 2: Drivers of change in immigration concern, 1-way FE mods, DKSEs")


# Table S3

form0 = diff(Imm_sup_z, lag=1) ~ plm::lag(Imm_sup_z, 1:2) + plm::lag(GDP_cp_grth, 1) + plm::lag(Unempl_WB, 1) 
form1 = update.formula(form0, . ~ . + ( diff(Imm_mus_pc, lag=1) + plm::lag(Imm_mus_pc, 1) ) * plm::lag(Noncit_stck_pc, 1))
form2 = update.formula(form1, . ~ . + plm::lag(FRight_seat_pr, 1) + plm::lag(BMI_i_ext, 1) + plm::lag(BMI_e_ext, 1) )
form0 = diff(Imm_sal_z, lag=1) ~ plm::lag(Imm_sal_z, 1:2) + plm::lag(GDP_cp_grth, 1) + plm::lag(Unempl_WB, 1) 
form3 = update.formula(form0, . ~ . + ( diff(Imm_mus_pc, lag=1) + plm::lag(Imm_mus_pc, 1) ) * plm::lag(Noncit_stck_pc, 1))
form4 = update.formula(form3, . ~ . + plm::lag(FRight_seat_pr, 1) + plm::lag(BMI_i_ext, 1) + plm::lag(BMI_e_ext, 1) )

mus.mods = lapply(list(form2, form4), plm, data=mig.plm, model="within", effect="indiv")
sapply(mus.mods, function(x) round(pwartest(x, order=1)[["p.value"]][[1]], 3) )
sapply(mus.mods, function(x) round(pcdtest(x, test="cd")[["p.value"]][[1]], 3) )
sapply(mus.mods, function(x) round(sigma(x), 3) )
vcm = lapply(mus.mods, function(x) plm::vcovSCC(x, cluster="time")) # cross-sectional heteroskedasticity and correlation
rse = lapply(vcm, function(x) sqrt(diag(x)))
pval = list()
for(i in 1:length(mus.mods)) {pval[[i]] = pt(abs(coef(mus.mods[[i]]) / rse[[i]]), df=500, lower.tail=FALSE) * 2}
screenreg(mus.mods, override.se=rse, override.pvalues=pval, digits=3, single.row=TRUE, leading.zero=FALSE, 
          stars=c(0.05))
vars = c("Immigration mood t-1", "Immigration mood t-2", "GDP growth per capita t-1", "Unemployment rate t-1",
         "Diff Muslim immigration inflows, t0", "Muslim immigration inflows, % of pop. t-1", 
         "Non-citizen stock, % of pop. t-1",  "Far right seat share t-1", "Immigrant integration policy index t-1", 
         "Immigrant entry policy index t-1", "Diff Muslim immigration inflows X non-citizen stock", 
         "Lag Muslim immigration inflows X non-citizen stock", "Immigration concern t-1", "Immigration concern t-2")
htmlreg(mus.mods, file="tableS3.htm", override.se=rse, override.pvalues=pval, digits=3, 
        single.row=TRUE, leading.zero=FALSE, stars=0.05, custom.coef.names=vars,
        caption="Table S3: Muslim immigration flows and changes in immigration mood and concern, 1-way FE mods, DKSEs")


# Table S2

form0 = diff(Imm_sup_z, lag=1) ~ plm::lag(Imm_sup_z, 1:2) + plm::lag(GDP_cp_grth, 1) + plm::lag(Unempl_WB, 1) + 
  plm::lag(FRight_seat_pr, 1) + plm::lag(BMI_i_ext, 1) + plm::lag(BMI_e_ext, 1)
form1 = update.formula(form0, . ~ . + ( diff(Net_migr_pc, lag=1) + plm::lag(Net_migr_pc, 1) ) * plm::lag(Noncit_stck_pc, 1))
form2 = update.formula(form0, . ~ . + ( diff(Immig_pc, lag=1) + plm::lag(Immig_pc, 1) ) * plm::lag(Noncit_stck_pc, 1))
form0 = diff(Imm_sal_z, lag=1) ~ plm::lag(Imm_sal_z, 1:2) + plm::lag(GDP_cp_grth, 1) + plm::lag(Unempl_WB, 1) + 
  plm::lag(FRight_seat_pr, 1) + plm::lag(BMI_i_ext, 1) + plm::lag(BMI_e_ext, 1)
form3 = update.formula(form0, . ~ . + ( diff(Net_migr_pc, lag=1) + plm::lag(Net_migr_pc, 1) ) * plm::lag(Noncit_stck_pc, 1))
form4 = update.formula(form0, . ~ . + ( diff(Immig_pc, lag=1) + plm::lag(Immig_pc, 1) ) * plm::lag(Noncit_stck_pc, 1))

twfe.mods = lapply(list(form1, form2, form3, form4), plm, data=mig.plm, model="within", effect="twoways")
sapply(twfe.mods, function(x) round(pwartest(x, order=1)[["p.value"]][[1]], 3) )
sapply(twfe.mods, function(x) round(pcdtest(x, test="cd")[["p.value"]][[1]], 3) )
sapply(twfe.mods, function(x) round(sigma(x), 3) )
vcm = lapply(twfe.mods, function(x) plm::vcovSCC(x, cluster="time")) # cross-sectional heteroskedasticity and correlation
rse = lapply(vcm, function(x) sqrt(diag(x)))
pval = list()
for(i in 1:length(twfe.mods)) {pval[[i]] = pt(abs(coef(twfe.mods[[i]]) / rse[[i]]), df=300, lower.tail=FALSE) * 2}
screenreg(twfe.mods, override.se=rse, override.pvalues=pval, digits=3, single.row=TRUE, leading.zero=FALSE, 
          stars=c(0.05))
vars = c("Immigration mood t-1", "Immigration mood t-2", "GDP growth per capita t-1", "Unemployment rate t-1",
         "Far right seat share t-1", "Immigrant integration policy index t-1", 
         "Immigrant entry policy index t-1", "Diff net migration, t0", "Net migration, t-1", 
         "Non-citizen stock, % of pop. t-1", "Diff net migration X non-citizen stock", 
         "Lag net migration X non-citizen stock", "Diff immigration inflows, t0", "Immigration inflows, t-1",
         "Diff immigration inflows X non-citizen stock", "Lag immigration inflows X non-citizen stock",
         "Immigration concern t-1", "Immigration concern t-2")
htmlreg(twfe.mods, file="tableS2.htm", override.se=rse, override.pvalues=pval, digits=3, 
        single.row=TRUE, leading.zero=FALSE, stars=0.05, custom.coef.names=vars,
        caption="Table S2: Drivers of change in immigration mood and concern, 2-way FE mods, DKSEs")


# Table S6: demographic models

imm.stck.mod.1 = plm(Noncit_stck_pc ~ plm::lag(Noncit_stck_pc, 1) + plm::lag(Net_migr_pc, 1), data=mig.plm, 
                     model="within", effect="indiv")
imm.stck.mod.2 = plm(Noncit_stck_pc ~ plm::lag(Noncit_stck_pc, 1) + plm::lag(Net_migr_pc, 1) + plm::lag(BMI_i_ext, 1) 
                     + plm::lag(Unempl_WB, 1), data=mig.plm, model="within", effect="indiv")
imm.stck.mod.3 = plm(Noncit_stck_pc ~ plm::lag(Noncit_stck_pc, 1) + plm::lag(Immig_pc, 1), data=mig.plm, model="within", 
                     effect="indiv")
imm.stck.mod.4 = plm(Noncit_stck_pc ~ plm::lag(Noncit_stck_pc, 1) + plm::lag(Immig_pc, 1) + plm::lag(BMI_i_ext, 1) 
                     + plm::lag(Unempl_WB, 1), data=mig.plm, model="within", effect="indiv")
dem.mods = list(imm.stck.mod.1, imm.stck.mod.2, imm.stck.mod.3, imm.stck.mod.4)

screenreg(dem.mods, single.row=TRUE, digits=3)

sapply(dem.mods, function(x) round(pwartest(x, order=1)[["p.value"]][[1]], 3) )
sapply(dem.mods, function(x) round(pcdtest(x, test="cd")[["p.value"]][[1]], 3) )
sapply(dem.mods, function(x) round(sigma(x), 3) )
vcm = lapply(dem.mods, function(x) plm::vcovSCC(x, cluster="time")) # cross-sectional heteroskedasticity and correlation
rse = lapply(vcm, function(x) sqrt(diag(x)))
pval = list()
for(i in 1:length(dem.mods)) {pval[[i]] = pt(abs(coef(dem.mods[[i]]) / rse[[i]]), df=500, lower.tail=FALSE) * 2}
screenreg(dem.mods, override.se=rse, override.pvalues=pval, digits=3, single.row=TRUE, leading.zero=FALSE, 
          stars=c(0.05))
vars = c("Non-citizen stock, % of pop. t-1", "Net migration, t-1", "Immigrant integration policy index t-1",
         "Unemployment rate t-1", "Immigration inflows t-1")
htmlreg(dem.mods, file="tableS6.htm", override.se=rse, override.pvalues=pval, digits=3, 
        single.row=TRUE, leading.zero=FALSE, stars=0.05, custom.coef.names=vars,
        caption="Table S6: Models of non-citizen stock accumulation, 1-way FE mods, DKSEs")


# Table S4, western europe only

form0 = diff(Imm_sup_z, lag=1) ~ plm::lag(Imm_sup_z, 1:2) + plm::lag(GDP_cp_grth, 1) + plm::lag(Unempl_WB, 1) + 
  plm::lag(FRight_seat_pr, 1) + plm::lag(BMI_i_ext, 1) + plm::lag(BMI_e_ext, 1)
form1 = update.formula(form0, . ~ . + ( diff(Net_migr_pc, lag=1) + plm::lag(Net_migr_pc, 1) ) * plm::lag(Noncit_stck_pc, 1))
form2 = update.formula(form0, . ~ . + ( diff(Immig_pc, lag=1) + plm::lag(Immig_pc, 1) ) * plm::lag(Noncit_stck_pc, 1))
form0 = diff(Imm_sal_z, lag=1) ~ plm::lag(Imm_sal_z, 1:2) + plm::lag(GDP_cp_grth, 1) + plm::lag(Unempl_WB, 1) + 
  plm::lag(FRight_seat_pr, 1) + plm::lag(BMI_i_ext, 1) + plm::lag(BMI_e_ext, 1)
form3 = update.formula(form0, . ~ . + ( diff(Net_migr_pc, lag=1) + plm::lag(Net_migr_pc, 1) ) * plm::lag(Noncit_stck_pc, 1))
form4 = update.formula(form0, . ~ . + ( diff(Immig_pc, lag=1) + plm::lag(Immig_pc, 1) ) * plm::lag(Noncit_stck_pc, 1))

we.mods = lapply(list(form1, form2, form3, form4), plm, data=mig.plm.we, model="within")
sapply(we.mods, function(x) round(pwartest(x, order=1)[["p.value"]][[1]], 3) )
sapply(we.mods, function(x) round(pcdtest(x, test="cd")[["p.value"]][[1]], 3) )
sapply(we.mods, function(x) round(sigma(x), 3) )
vcm = lapply(we.mods, function(x) plm::vcovSCC(x, cluster="time")) # cross-sectional heteroskedasticity and correlation
rse = lapply(vcm, function(x) sqrt(diag(x)))
pval = list()
for(i in 1:length(we.mods)) {pval[[i]] = pt(abs(coef(we.mods[[i]]) / rse[[i]]), df=300, lower.tail=FALSE) * 2}
screenreg(we.mods, override.se=rse, override.pvalues=pval, digits=3, single.row=TRUE, leading.zero=FALSE, 
          stars=c(0.05))
vars = c("Immigration mood t-1", "Immigration mood t-2", "GDP growth per capita t-1", "Unemployment rate t-1",
         "Far right seat share t-1", "Immigrant integration policy index t-1", 
         "Immigrant entry policy index t-1", "Diff net migration, t0", "Net migration, t-1", 
         "Non-citizen stock, % of pop. t-1", "Diff net migration X non-citizen stock", 
         "Lag net migration X non-citizen stock", "Diff immigration inflows, t0", "Immigration inflows, t-1",
         "Diff immigration inflows X non-citizen stock", "Lag immigration inflows X non-citizen stock",
         "Immigration concern t-1", "Immigration concern t-2")
htmlreg(we.mods, file="tableS4.htm", override.se=rse, override.pvalues=pval, digits=3, 
        single.row=TRUE, leading.zero=FALSE, stars=0.05, custom.coef.names=vars,
        caption="Table S4: Drivers of change in immigration mood and concern, Western Europe only, DKSEs")


# Table S5: citizen survey samples only

form0 = diff(Imm_sup_cit_z, lag=1) ~ plm::lag(Imm_sup_cit_z, 1:2) + plm::lag(GDP_cp_grth, 1) + plm::lag(Unempl_WB, 1) + 
  plm::lag(FRight_seat_pr, 1) + plm::lag(BMI_i_ext, 1) + plm::lag(BMI_e_ext, 1)
form1 = update.formula(form0, . ~ . + ( diff(Net_migr_pc, lag=1) + plm::lag(Net_migr_pc, 1) ) * plm::lag(Noncit_stck_pc, 1))
form2 = update.formula(form0, . ~ . + ( diff(Immig_pc, lag=1) + plm::lag(Immig_pc, 1) ) * plm::lag(Noncit_stck_pc, 1))

cit.mods = lapply(list(form1, form2), plm, data=mig.plm, model="within")
sapply(cit.mods, function(x) round(pwartest(x, order=1)[["p.value"]][[1]], 3) )
sapply(cit.mods, function(x) round(pcdtest(x, test="cd")[["p.value"]][[1]], 3) )
sapply(cit.mods, function(x) round(sigma(x), 3) )
vcm = lapply(cit.mods, function(x) plm::vcovSCC(x, cluster="time")) # cross-sectional heteroskedasticity and correlation
rse = lapply(vcm, function(x) sqrt(diag(x)))
pval = list()
for(i in 1:length(cit.mods)) {pval[[i]] = pt(abs(coef(cit.mods[[i]]) / rse[[i]]), df=300, lower.tail=FALSE) * 2}
screenreg(cit.mods, override.se=rse, override.pvalues=pval, digits=3, single.row=TRUE, leading.zero=FALSE, 
          stars=c(0.05))
vars = c("Immigration mood t-1", "Immigration mood t-2", "GDP growth per capita t-1", "Unemployment rate t-1",
         "Far right seat share t-1", "Immigrant integration policy index t-1", 
         "Immigrant entry policy index t-1", "Diff net migration, t0", "Net migration, t-1", 
         "Non-citizen stock, % of pop. t-1", "Diff net migration X non-citizen stock", 
         "Lag net migration X non-citizen stock", "Diff immigration inflows, t0", "Immigration inflows, t-1",
         "Diff immigration inflows X non-citizen stock", "Lag immigration inflows X non-citizen stock")
htmlreg(cit.mods, file="tableS5.htm", override.se=rse, override.pvalues=pval, digits=3, 
        single.row=TRUE, leading.zero=FALSE, stars=0.05, custom.coef.names=vars,
        caption="Table S5: Drivers of change in immigration mood, citizen samples, DKSEs")



### Marginal effects

# Figure S3a

marg.min = min(migrat$Noncit_stck_within, na.rm=TRUE)
marg.max = max(migrat$Noncit_stck_within, na.rm=TRUE)

marg.mod = perc.mods[[2]] 
marg.vcov = plm::vcovSCC(marg.mod, cluster="time")
sim.stck.data = seq(from=marg.min, to=marg.max, length.out=100)
me.pe = coef(marg.mod)["plm::lag(Immig_pc, 1)"] + 
        (coef(marg.mod)["plm::lag(Immig_pc, 1):plm::lag(Noncit_stck_pc, 1)"] * sim.stck.data)
me.se = sqrt(marg.vcov["plm::lag(Immig_pc, 1)", "plm::lag(Immig_pc, 1)"] + 
               (sim.stck.data^2 * marg.vcov["plm::lag(Immig_pc, 1):plm::lag(Noncit_stck_pc, 1)", 
                                            "plm::lag(Immig_pc, 1):plm::lag(Noncit_stck_pc, 1)"]) + 
               (2 * sim.stck.data * marg.vcov["plm::lag(Immig_pc, 1)", 
                                              "plm::lag(Immig_pc, 1):plm::lag(Noncit_stck_pc, 1)"]))

dens.stck = density(migrat$Noncit_stck_within, na.rm=TRUE)
dens.stck$y = (dens.stck$y / 3) - 0.3
pdf("figS3a.pdf", height=4, width=4)
par(mfrow=c(1,1), mar=c(2.8,3,0.5,0.2), tcl=-0.3, las=0, cex=0.9, mgp=c(1,0.3,0))
plot(x=sim.stck.data, y=me.pe, type="n", ylim=c(-0.3, 0.3), xaxs="i", yaxs="i", xlab="", ylab="", 
     cex.axis=0.85)
grid(col=rgb(0,0,0,0.5), lwd=0.75)
abline(h=0, lty=1, col=rgb(0,0,0,0.7))
lines(x=sim.stck.data, y=me.pe, lty=1, lwd=2, col=rgb(0,0.4,0.4,1))
polygon(x=c(sim.stck.data, sim.stck.data[100:1]), y=c((me.pe - 1.96*me.se), (me.pe + 1.96*me.se)[100:1]), 
        col=rgb(0,0.4,0.4,0.4), border=NA)
polygon(x=c(dens.stck$x, sim.stck.data[100:1]), y=c(dens.stck$y, rep(-0.3, 100)), col=rgb(0,0,0,0.4), 
        border=rgb(0,0,0,0.8))
mtext(side=1, line=1.3, text="Non-citizens as % of population, within-country", cex=0.9)
mtext(side=2, line=1.5, text="Marginal effects of immigration inflows", cex=0.9)
text(x=-6.2, y=0.28, label="Immigration mood", cex=1, adj=0)
dev.off()


# Figure S3b

marg.mod = perc.mods[[1]] 
marg.vcov = plm::vcovSCC(marg.mod, cluster="time")
sim.stck.data = seq(from=marg.min, to=marg.max, length.out=100)
me.pe = coef(marg.mod)["plm::lag(Net_migr_pc, 1)"] + 
        (coef(marg.mod)["plm::lag(Net_migr_pc, 1):plm::lag(Noncit_stck_pc, 1)"] * sim.stck.data)
me.se = sqrt(marg.vcov["plm::lag(Net_migr_pc, 1)", "plm::lag(Net_migr_pc, 1)"] + 
               (sim.stck.data^2 * marg.vcov["plm::lag(Net_migr_pc, 1):plm::lag(Noncit_stck_pc, 1)", 
                                            "plm::lag(Net_migr_pc, 1):plm::lag(Noncit_stck_pc, 1)"]) + 
               (2 * sim.stck.data * marg.vcov["plm::lag(Net_migr_pc, 1)", 
                                              "plm::lag(Net_migr_pc, 1):plm::lag(Noncit_stck_pc, 1)"]))

dens.stck = density(migrat$Noncit_stck_within, na.rm=TRUE)
dens.stck$y = (dens.stck$y / 3) - 0.3
pdf("figS3b.pdf", height=4, width=4)
par(mfrow=c(1,1), mar=c(2.8,3,0.5,0.2), tcl=-0.3, las=0, cex=0.9, mgp=c(1,0.3,0))
plot(x=sim.stck.data, y=me.pe, type="n", ylim=c(-0.3, 0.3), xaxs="i", yaxs="i", xlab="", ylab="", 
     cex.axis=0.85)
grid(col=rgb(0,0,0,0.5), lwd=0.75)
abline(h=0, lty=1, col=rgb(0,0,0,0.7))
lines(x=sim.stck.data, y=me.pe, lty=1, lwd=2, col=rgb(0,0.4,0.4,1))
polygon(x=c(sim.stck.data, sim.stck.data[100:1]), y=c((me.pe - 1.96*me.se), (me.pe + 1.96*me.se)[100:1]), 
        col=rgb(0,0.4,0.4,0.4), border=NA)
polygon(x=c(dens.stck$x, sim.stck.data[100:1]), y=c(dens.stck$y, rep(-0.3, 100)), col=rgb(0,0,0,0.4), 
        border=rgb(0,0,0,0.8))
mtext(side=1, line=1.3, text="Non-citizens as % of population, within-country", cex=0.9)
mtext(side=2, line=1.5, text="Marginal effects of net migration", cex=0.9)
text(x=-6.2, y=0.28, label="Immigration mood", cex=1, adj=0)
dev.off()




#### Simulated Effects Plots 
#### These project future immigrant stock based on current immigration flows

## Figure 2a

imm.stck.mod.1 = plm(Noncit_stck_pc ~ plm::lag(Noncit_stck_pc, 1) + plm::lag(Net_migr_pc, 1), data=mig.plm, model="within", effect="indiv")

mod = perc.mods[[1]]
dem.mod = imm.stck.mod.1
round(coef(mod), 2)
n.coef =  length(coef(mod))
n.yrs = 30
n.burn = 0

summary(migrat$Net_migr_within)
sd(migrat$Net_migr_within, na.rm=TRUE)
summary(migrat$Noncit_stck_within)
sd(migrat$Noncit_stck_within, na.rm=TRUE)

chg.nm = 2 * sd(migrat$Net_migr_within, na.rm=TRUE)
pred.data = data.frame(ImmSup_m1=0, ImmSup_m2=0, 
                       GDPgr=0, Unemp=0, 
                       diff_NM=0,
                       NM=c(rep(0, n.burn), rep(chg.nm, n.yrs-n.burn)), 
                       Noncit_stck_pc=-5,
                       FRseat=0, BMIi=0, BMIe=0,
                       DiffFlowsxStock=NA,
                       FlowsxStock=NA,
                       ChgSup=NA, ImmSup=NA
                       )
n.sims = 1e3
sim.data = array(as.matrix(pred.data), dim=c(dim(pred.data), n.sims))
mod.sims = mvrnorm(n=n.sims, mu=coef(mod), Sigma=plm::vcovSCC(mod, cluster="time"), empirical=FALSE)
dem.mod.sims = mvrnorm(n=n.sims, mu=coef(dem.mod), Sigma=plm::vcovSCC(dem.mod, cluster="time"), empirical=FALSE)
dem.hat = matrix(NA, ncol=n.yrs, nrow=n.sims)
yhat = matrix(NA, ncol=n.yrs, nrow=n.sims)
sup.hat = matrix(NA, ncol=n.yrs, nrow=n.sims)

for(i in 1:(n.yrs-1)) {
  for(k in 1:n.sims) {
    dem.hat[k, i] = dem.mod.sims[k,] %*% as.matrix(sim.data[i, c(7,6), k])
    sim.data[i+1, match("Noncit_stck_pc", names(pred.data)), k] = dem.hat[k, i] + rnorm(1, 0, sigma(dem.mod))
    sim.data[i, match("FlowsxStock", names(pred.data)), k] = sim.data[i, match("NM", names(pred.data)), k] * 
      sim.data[i, match("Noncit_stck_pc", names(pred.data)), k]
    sim.data[i+1, match("diff_NM", names(pred.data)), k] = sim.data[i+1, match("NM", names(pred.data)), k] - 
      sim.data[i, match("NM", names(pred.data)), k]
    sim.data[i, match("DiffFlowsxStock", names(pred.data)), k] = sim.data[i, match("diff_NM", names(pred.data)), k] * 
      sim.data[i, match("Noncit_stck_pc", names(pred.data)), k]
    sim.data[i, match("ChgSup", names(pred.data)), k] = sum(sim.data[i, 1:n.coef, k] * coef(mod))
    sim.data[i, match("ImmSup", names(pred.data)), k] = sim.data[i, match("ImmSup_m1", names(pred.data)), k] + 
      sim.data[i, match("ChgSup", names(pred.data)), k]
    sim.data[i+1, match("ImmSup_m1", names(pred.data)), k] = sim.data[i, match("ImmSup", names(pred.data)), k]
    sim.data[i+1, match("ImmSup_m2", names(pred.data)), k] = sim.data[i, match("ImmSup_m1", names(pred.data)), k]    
    sim.data[i, match("ChgSup", names(pred.data)), k] = mod.sims[k,] %*% (as.matrix(sim.data[i, 1:n.coef, k]))
    sup.hat[k, i] = sim.data[i, match("ImmSup_m1", names(pred.data)), k] + sim.data[i,match("ChgSup", names(pred.data)), k]
  }
}

pe = apply(sup.hat, 2, mean, na.rm=TRUE)
pe.sd = apply(sup.hat, 2, sd, na.rm=TRUE)
u95 = pe + 1.96 * pe.sd
l95 = pe - 1.96 * pe.sd

yr.st = 1
yr1.off = pe[yr.st]
pe = pe-yr1.off
u95 = u95-yr1.off
l95 = l95-yr1.off

pdf("fig2a.pdf", height=4, width=4)
par(mfrow=c(1,1), mar=c(2,3,2,.5), tcl=-0.25, cex=0.9)
plot(x=yr.st:(n.yrs-1), y=pe[yr.st:(n.yrs-1)], type="n", ylim=c(-1, 1), 
     xaxs="i", yaxs="i", xaxt="n", yaxt="n", xlab="", ylab="")
abline(v=seq(10, 30, by=10), lty=3, col=rgb(0,0,0,.4))
abline(h=c(-0.5, 0, 0.5), lty=3, col=rgb(0,0,0,.4))
abline(v=0, lty=2, col=rgb(0,0,0,.75))
axis(side=1, at=seq(0, 30, by=10), labels=seq(0, 30, by=10), cex.axis=0.8, las=1, mgp=c(0.7, 0.2, 0))
axis(side=2, at=c(-1,-0.5,0,0.5,1), cex.axis=0.8, las=0, mgp=c(1, 0.4, 0))
lines(x=yr.st:(n.yrs-1), y=pe[yr.st:(n.yrs-1)], col=rgb(0,.4,.4,1), lwd=2)
polygon(x=c(yr.st:(n.yrs-1), (n.yrs-1):yr.st), col=rgb(0,.4,.4,.5), border=NA, 
        y=c(u95[yr.st:(n.yrs-1)], l95[(n.yrs-1):yr.st]))
mtext(side=1, line=1.1, "Years", cex=0.9)
mtext(side=2, line=1.6, "Immigration mood", cex=0.9)
text(x=0.5, y=0.8, "Net migration increases by 2 SD in year 0", cex=0.85, pos=4)
mtext(text="Existing immigrant stock is low", side=3, line=0.4, cex=0.9, adj=0)
dev.off()


## Figure 3a

mod = sal.mods[[1]]
dem.mod = imm.stck.mod.1
round(coef(mod), 2)
n.coef =  length(coef(mod))
n.yrs = 30
n.burn = 0

summary(migrat$Net_migr_within)
sd(migrat$Net_migr_within, na.rm=TRUE)
summary(migrat$Noncit_stck_within)
sd(migrat$Noncit_stck_within, na.rm=TRUE)

chg.nm = 2 * sd(migrat$Net_migr_within, na.rm=TRUE)
pred.data = data.frame(ImmSal_m1=0, ImmSal_m2=0, 
                       GDPgr=0, Unemp=0, 
                       diff_NM=0,
                       NM=c(rep(0, n.burn), rep(chg.nm, n.yrs-n.burn)), 
                       Noncit_stck_pc=-5,
                       FRseat=0, BMIi=0, BMIe=0,
                       DiffFlowsxStock=NA,
                       FlowsxStock=NA,
                       ChgSal=NA, ImmSal=NA
                       )
n.sims = 1e3
sim.data = array(as.matrix(pred.data), dim=c(dim(pred.data), n.sims))
mod.sims = mvrnorm(n=n.sims, mu=coef(mod), Sigma=plm::vcovSCC(mod, cluster="time"), empirical=FALSE)
dem.mod.sims = mvrnorm(n=n.sims, mu=coef(dem.mod), Sigma=plm::vcovSCC(dem.mod, cluster="time"), empirical=FALSE)
dem.hat = matrix(NA, ncol=n.yrs, nrow=n.sims)
yhat = matrix(NA, ncol=n.yrs, nrow=n.sims)
sup.hat = matrix(NA, ncol=n.yrs, nrow=n.sims)

for(i in 1:(n.yrs-1)) {
  for(k in 1:n.sims) {
    dem.hat[k, i] = dem.mod.sims[k,] %*% as.matrix(sim.data[i, c(7,6), k])
    sim.data[i+1, match("Noncit_stck_pc", names(pred.data)), k] = dem.hat[k, i] + rnorm(1, 0, sigma(dem.mod))
    sim.data[i, match("FlowsxStock", names(pred.data)), k] = sim.data[i, match("NM", names(pred.data)), k] * 
      sim.data[i, match("Noncit_stck_pc", names(pred.data)), k]
    sim.data[i+1, match("diff_NM", names(pred.data)), k] = sim.data[i+1, match("NM", names(pred.data)), k] - 
      sim.data[i, match("NM", names(pred.data)), k]
    sim.data[i, match("DiffFlowsxStock", names(pred.data)), k] = sim.data[i, match("diff_NM", names(pred.data)), k] * 
      sim.data[i, match("Noncit_stck_pc", names(pred.data)), k]
    sim.data[i, match("ChgSal", names(pred.data)), k] = sum(sim.data[i, 1:n.coef, k] * coef(mod))
    sim.data[i, match("ImmSal", names(pred.data)), k] = sim.data[i, match("ImmSal_m1", names(pred.data)), k] + 
      sim.data[i, match("ChgSal", names(pred.data)), k]
    sim.data[i+1, match("ImmSal_m1", names(pred.data)), k] = sim.data[i, match("ImmSal", names(pred.data)), k]
    sim.data[i+1, match("ImmSal_m2", names(pred.data)), k] = sim.data[i, match("ImmSal_m1", names(pred.data)), k]    
    sim.data[i, match("ChgSal", names(pred.data)), k] = mod.sims[k,] %*% (as.matrix(sim.data[i, 1:n.coef, k]))
    sup.hat[k, i] = sim.data[i, match("ImmSal_m1", names(pred.data)), k] + sim.data[i,match("ChgSal", names(pred.data)), k]
  }
}

pe = apply(sup.hat, 2, mean, na.rm=TRUE)
pe.sd = apply(sup.hat, 2, sd, na.rm=TRUE)
u95 = pe + 1.96 * pe.sd
l95 = pe - 1.96 * pe.sd

yr.st = 1
yr1.off = pe[yr.st]
pe = pe-yr1.off
u95 = u95-yr1.off
l95 = l95-yr1.off

pdf("fig3a.pdf", height=4, width=4)
par(mfrow=c(1,1), mar=c(2,3,2,.5), tcl=-0.25, cex=0.9)
plot(x=yr.st:(n.yrs-1), y=pe[yr.st:(n.yrs-1)], type="n", ylim=c(-1, 1), 
     xaxs="i", yaxs="i", xaxt="n", yaxt="n", xlab="", ylab="")
abline(v=seq(10, 30, by=10), lty=3, col=rgb(0,0,0,.4))
abline(h=c(-0.5, 0, 0.5), lty=3, col=rgb(0,0,0,.4))
abline(v=0, lty=2, col=rgb(0,0,0,.75))
axis(side=1, at=seq(0, 30, by=10), labels=seq(0, 30, by=10), cex.axis=0.8, las=1, mgp=c(0.7, 0.2, 0))
axis(side=2, at=c(-1,-0.5,0,0.5,1), cex.axis=0.8, las=0, mgp=c(1, 0.4, 0))
lines(x=yr.st:(n.yrs-1), y=pe[yr.st:(n.yrs-1)], col=rgb(0,.4,.4,1), lwd=2)
polygon(x=c(yr.st:(n.yrs-1), (n.yrs-1):yr.st), col=rgb(0,.4,.4,.5), border=NA, 
        y=c(u95[yr.st:(n.yrs-1)], l95[(n.yrs-1):yr.st]))
mtext(side=1, line=1.1, "Years", cex=0.9)
mtext(side=2, line=1.6, "Immigration concern", cex=0.9)
text(x=0.5, y=0.8, "Net migration increases by 2 SD in year 0", cex=0.85, pos=4)
mtext(text="Existing immigrant stock is low", side=3, line=0.4, cex=0.9, adj=0)
dev.off()


## Figure 2b

imm.stck.mod.1 = plm(Noncit_stck_pc ~ plm::lag(Noncit_stck_pc, 1) + plm::lag(Immig_pc, 1), data=mig.plm, model="within", effect="indiv")

mod = perc.mods[[2]]
dem.mod = imm.stck.mod.1
round(coef(mod), 2)
n.coef =  length(coef(mod))
n.yrs = 30
n.burn = 0

summary(migrat$Immig_pc_within)
sd(migrat$Immig_pc_within, na.rm=TRUE)
summary(migrat$Noncit_stck_within)
sd(migrat$Noncit_stck_within, na.rm=TRUE)

chg.immin = 2 * sd(migrat$Immig_pc_within, na.rm=TRUE)
pred.data = data.frame(ImmSup_m1=0, ImmSup_m2=0, 
                       GDPgr=0, Unemp=0, 
                       diff_ImmIn=0,
                       ImmIn=c(rep(0, n.burn), rep(chg.immin, n.yrs-n.burn)), 
                       Noncit_stck_pc=-5,
                       FRseat=0, BMIi=0, BMIe=0,
                       DiffFlowsxStock=NA,
                       FlowsxStock=NA,
                       ChgSup=NA, ImmSup=NA
)
n.sims = 1e3
sim.data = array(as.matrix(pred.data), dim=c(dim(pred.data), n.sims))
mod.sims = mvrnorm(n=n.sims, mu=coef(mod), Sigma=plm::vcovSCC(mod, cluster="time"), empirical=FALSE)
dem.mod.sims = mvrnorm(n=n.sims, mu=coef(dem.mod), Sigma=plm::vcovSCC(dem.mod, cluster="time"), empirical=FALSE)
dem.hat = matrix(NA, ncol=n.yrs, nrow=n.sims)
yhat = matrix(NA, ncol=n.yrs, nrow=n.sims)
sup.hat = matrix(NA, ncol=n.yrs, nrow=n.sims)

for(i in 1:(n.yrs-1)) {
  for(k in 1:n.sims) {
    dem.hat[k, i] = dem.mod.sims[k,] %*% as.matrix(sim.data[i, c(7,6), k])
    sim.data[i+1, match("Noncit_stck_pc", names(pred.data)), k] = dem.hat[k, i] + rnorm(1, 0, sigma(dem.mod))
    sim.data[i, match("FlowsxStock", names(pred.data)), k] = sim.data[i, match("ImmIn", names(pred.data)), k] * 
      sim.data[i, match("Noncit_stck_pc", names(pred.data)), k]
    sim.data[i+1, match("diff_ImmIn", names(pred.data)), k] = sim.data[i+1, match("ImmIn", names(pred.data)), k] - 
      sim.data[i, match("ImmIn", names(pred.data)), k]
    sim.data[i, match("DiffFlowsxStock", names(pred.data)), k] = sim.data[i, match("diff_ImmIn", names(pred.data)), k] * 
      sim.data[i, match("Noncit_stck_pc", names(pred.data)), k]
    sim.data[i, match("ChgSup", names(pred.data)), k] = sum(sim.data[i, 1:n.coef, k] * coef(mod))
    sim.data[i, match("ImmSup", names(pred.data)), k] = sim.data[i, match("ImmSup_m1", names(pred.data)), k] + 
      sim.data[i, match("ChgSup", names(pred.data)), k]
    sim.data[i+1, match("ImmSup_m1", names(pred.data)), k] = sim.data[i, match("ImmSup", names(pred.data)), k]
    sim.data[i+1, match("ImmSup_m2", names(pred.data)), k] = sim.data[i, match("ImmSup_m1", names(pred.data)), k]    
    sim.data[i, match("ChgSup", names(pred.data)), k] = mod.sims[k,] %*% (as.matrix(sim.data[i, 1:n.coef, k]))
    sup.hat[k, i] = sim.data[i, match("ImmSup_m1", names(pred.data)), k] + sim.data[i,match("ChgSup", names(pred.data)), k]
  }
}

pe = apply(sup.hat, 2, mean, na.rm=TRUE)
pe.sd = apply(sup.hat, 2, sd, na.rm=TRUE)
u95 = pe + 1.96 * pe.sd
l95 = pe - 1.96 * pe.sd

yr.st = 1
yr1.off = pe[yr.st]
pe = pe-yr1.off
u95 = u95-yr1.off
l95 = l95-yr1.off

pdf("fig2b.pdf", height=4, width=4)
par(mfrow=c(1,1), mar=c(2,3,2,.5), tcl=-0.25, cex=0.9)
plot(x=yr.st:(n.yrs-1), y=pe[yr.st:(n.yrs-1)], type="n", ylim=c(-1, 1), 
     xaxs="i", yaxs="i", xaxt="n", yaxt="n", xlab="", ylab="")
abline(v=seq(10, 30, by=10), lty=3, col=rgb(0,0,0,.4))
abline(h=c(-0.5, 0, 0.5), lty=3, col=rgb(0,0,0,.4))
abline(v=0, lty=2, col=rgb(0,0,0,.75))
axis(side=1, at=seq(0, 30, by=10), labels=seq(0, 30, by=10), cex.axis=0.8, las=1, mgp=c(0.7, 0.2, 0))
axis(side=2, at=c(-1,-0.5,0,0.5,1), cex.axis=0.8, las=0, mgp=c(1, 0.4, 0))
lines(x=yr.st:(n.yrs-1), y=pe[yr.st:(n.yrs-1)], col=rgb(0,.4,.4,1), lwd=2)
polygon(x=c(yr.st:(n.yrs-1), (n.yrs-1):yr.st), col=rgb(0,.4,.4,.5), border=NA, 
        y=c(u95[yr.st:(n.yrs-1)], l95[(n.yrs-1):yr.st]))
mtext(side=1, line=1.1, "Years", cex=0.9)
mtext(side=2, line=1.6, "Immigration mood", cex=0.9)
text(x=0.5, y=0.8, "Immigration inflows increase by 2 SD in year 0", cex=0.85, pos=4)
mtext(text="Existing immigrant stock is low", side=3, line=0.4, cex=0.9, adj=0)
dev.off()


## Figure 3b

mod = sal.mods[[2]]
dem.mod = imm.stck.mod.1
round(coef(mod), 2)
n.coef =  length(coef(mod))
n.yrs = 30
n.burn = 0

summary(migrat$Immig_pc_within)
sd(migrat$Immig_pc_within, na.rm=TRUE)
summary(migrat$Noncit_stck_within)
sd(migrat$Noncit_stck_within, na.rm=TRUE)

chg.immin = 2 * sd(migrat$Immig_pc_within, na.rm=TRUE)
pred.data = data.frame(ImmSal_m1=0, ImmSal_m2=0, 
                       GDPgr=0, Unemp=0, 
                       diff_ImmIn=0,
                       ImmIn=c(rep(0, n.burn), rep(chg.immin, n.yrs-n.burn)), 
                       Noncit_stck_pc=-5,
                       FRseat=0, BMIi=0, BMIe=0,
                       DiffFlowsxStock=NA,
                       FlowsxStock=NA,
                       ChgSal=NA, ImmSal=NA
)
n.sims = 1e3
sim.data = array(as.matrix(pred.data), dim=c(dim(pred.data), n.sims))
mod.sims = mvrnorm(n=n.sims, mu=coef(mod), Sigma=plm::vcovSCC(mod, cluster="time"), empirical=FALSE)
dem.mod.sims = mvrnorm(n=n.sims, mu=coef(dem.mod), Sigma=plm::vcovSCC(dem.mod, cluster="time"), empirical=FALSE)
dem.hat = matrix(NA, ncol=n.yrs, nrow=n.sims)
yhat = matrix(NA, ncol=n.yrs, nrow=n.sims)
sup.hat = matrix(NA, ncol=n.yrs, nrow=n.sims)

for(i in 1:(n.yrs-1)) {
  for(k in 1:n.sims) {
    dem.hat[k, i] = dem.mod.sims[k,] %*% as.matrix(sim.data[i, c(7,6), k])
    sim.data[i+1, match("Noncit_stck_pc", names(pred.data)), k] = dem.hat[k, i] + rnorm(1, 0, sigma(dem.mod))
    sim.data[i, match("FlowsxStock", names(pred.data)), k] = sim.data[i, match("ImmIn", names(pred.data)), k] * 
      sim.data[i, match("Noncit_stck_pc", names(pred.data)), k]
    sim.data[i+1, match("diff_ImmIn", names(pred.data)), k] = sim.data[i+1, match("ImmIn", names(pred.data)), k] - 
      sim.data[i, match("ImmIn", names(pred.data)), k]
    sim.data[i, match("DiffFlowsxStock", names(pred.data)), k] = sim.data[i, match("diff_ImmIn", names(pred.data)), k] * 
      sim.data[i, match("Noncit_stck_pc", names(pred.data)), k]
    sim.data[i, match("ChgSal", names(pred.data)), k] = sum(sim.data[i, 1:n.coef, k] * coef(mod))
    sim.data[i, match("ImmSal", names(pred.data)), k] = sim.data[i, match("ImmSal_m1", names(pred.data)), k] + 
      sim.data[i, match("ChgSal", names(pred.data)), k]
    sim.data[i+1, match("ImmSal_m1", names(pred.data)), k] = sim.data[i, match("ImmSal", names(pred.data)), k]
    sim.data[i+1, match("ImmSal_m2", names(pred.data)), k] = sim.data[i, match("ImmSal_m1", names(pred.data)), k]    
    sim.data[i, match("ChgSal", names(pred.data)), k] = mod.sims[k,] %*% (as.matrix(sim.data[i, 1:n.coef, k]))
    sup.hat[k, i] = sim.data[i, match("ImmSal_m1", names(pred.data)), k] + sim.data[i,match("ChgSal", names(pred.data)), k]
  }
}

pe = apply(sup.hat, 2, mean, na.rm=TRUE)
pe.sd = apply(sup.hat, 2, sd, na.rm=TRUE)
u95 = pe + 1.96 * pe.sd
l95 = pe - 1.96 * pe.sd

yr.st = 1
yr1.off = pe[yr.st]
pe = pe-yr1.off
u95 = u95-yr1.off
l95 = l95-yr1.off

pdf("fig3b.pdf", height=4, width=4)
par(mfrow=c(1,1), mar=c(2,3,2,.5), tcl=-0.25, cex=0.9)
plot(x=yr.st:(n.yrs-1), y=pe[yr.st:(n.yrs-1)], type="n", ylim=c(-1, 1), 
     xaxs="i", yaxs="i", xaxt="n", yaxt="n", xlab="", ylab="")
abline(v=seq(10, 30, by=10), lty=3, col=rgb(0,0,0,.4))
abline(h=c(-0.5, 0, 0.5), lty=3, col=rgb(0,0,0,.4))
abline(v=0, lty=2, col=rgb(0,0,0,.75))
axis(side=1, at=seq(0, 30, by=10), labels=seq(0, 30, by=10), cex.axis=0.8, las=1, mgp=c(0.7, 0.2, 0))
axis(side=2, at=c(-1,-0.5,0,0.5,1), cex.axis=0.8, las=0, mgp=c(1, 0.4, 0))
lines(x=yr.st:(n.yrs-1), y=pe[yr.st:(n.yrs-1)], col=rgb(0,.4,.4,1), lwd=2)
polygon(x=c(yr.st:(n.yrs-1), (n.yrs-1):yr.st), col=rgb(0,.4,.4,.5), border=NA, 
        y=c(u95[yr.st:(n.yrs-1)], l95[(n.yrs-1):yr.st]))
mtext(side=1, line=1.1, "Years", cex=0.9)
mtext(side=2, line=1.6, "Immigration concern", cex=0.9)
text(x=0.5, y=0.8, "Immigration inflows increase by 2 SD in year 0", cex=0.85, pos=4)
mtext(text="Existing immigrant stock is low", side=3, line=0.4, cex=0.9, adj=0)
dev.off()


## Figure 2c

imm.stck.mod.1 = plm(Noncit_stck_pc ~ plm::lag(Noncit_stck_pc, 1) + plm::lag(Net_migr_pc, 1), data=mig.plm, model="within", effect="indiv")
mod = perc.mods[[1]]
dem.mod = imm.stck.mod.1
round(coef(mod), 2)
n.coef =  length(coef(mod))
n.yrs = 30
n.burn = 0

chg.nm = 2 * sd(migrat$Net_migr_within, na.rm=TRUE)
pred.data = data.frame(ImmSup_m1=0, ImmSup_m2=0, 
                       GDPgr=0, Unemp=0, 
                       diff_NM=0,
                       NM=c(rep(0, n.burn), rep(chg.nm, n.yrs-n.burn)), 
                       Noncit_stck_pc=0,
                       FRseat=0, BMIi=0, BMIe=0,
                       DiffFlowsxStock=NA,
                       FlowsxStock=NA,
                       ChgSup=NA, ImmSup=NA
)
n.sims = 1e3
sim.data = array(as.matrix(pred.data), dim=c(dim(pred.data), n.sims))
mod.sims = mvrnorm(n=n.sims, mu=coef(mod), Sigma=plm::vcovSCC(mod, cluster="time"), empirical=FALSE)
dem.mod.sims = mvrnorm(n=n.sims, mu=coef(dem.mod), Sigma=plm::vcovSCC(dem.mod, cluster="time"), empirical=FALSE)
dem.hat = matrix(NA, ncol=n.yrs, nrow=n.sims)
yhat = matrix(NA, ncol=n.yrs, nrow=n.sims)
sup.hat = matrix(NA, ncol=n.yrs, nrow=n.sims)

for(i in 1:(n.yrs-1)) {
  for(k in 1:n.sims) {
    dem.hat[k, i] = dem.mod.sims[k,] %*% as.matrix(sim.data[i, c(7,6), k])
    sim.data[i+1, match("Noncit_stck_pc", names(pred.data)), k] = dem.hat[k, i] + rnorm(1, 0, sigma(dem.mod))
    sim.data[i, match("FlowsxStock", names(pred.data)), k] = sim.data[i, match("NM", names(pred.data)), k] * 
      sim.data[i, match("Noncit_stck_pc", names(pred.data)), k]
    sim.data[i+1, match("diff_NM", names(pred.data)), k] = sim.data[i+1, match("NM", names(pred.data)), k] - 
      sim.data[i, match("NM", names(pred.data)), k]
    sim.data[i, match("DiffFlowsxStock", names(pred.data)), k] = sim.data[i, match("diff_NM", names(pred.data)), k] * 
      sim.data[i, match("Noncit_stck_pc", names(pred.data)), k]
    sim.data[i, match("ChgSup", names(pred.data)), k] = sum(sim.data[i, 1:n.coef, k] * coef(mod))
    sim.data[i, match("ImmSup", names(pred.data)), k] = sim.data[i, match("ImmSup_m1", names(pred.data)), k] + 
      sim.data[i, match("ChgSup", names(pred.data)), k]
    sim.data[i+1, match("ImmSup_m1", names(pred.data)), k] = sim.data[i, match("ImmSup", names(pred.data)), k]
    sim.data[i+1, match("ImmSup_m2", names(pred.data)), k] = sim.data[i, match("ImmSup_m1", names(pred.data)), k]    
    sim.data[i, match("ChgSup", names(pred.data)), k] = mod.sims[k,] %*% (as.matrix(sim.data[i, 1:n.coef, k]))
    sup.hat[k, i] = sim.data[i, match("ImmSup_m1", names(pred.data)), k] + sim.data[i,match("ChgSup", names(pred.data)), k]
  }
}

pe = apply(sup.hat, 2, mean, na.rm=TRUE)
pe.sd = apply(sup.hat, 2, sd, na.rm=TRUE)
u95 = pe + 1.96 * pe.sd
l95 = pe - 1.96 * pe.sd
yr.st = 1
yr1.off = pe[yr.st]
pe = pe-yr1.off
u95 = u95-yr1.off
l95 = l95-yr1.off

pdf("fig2c.pdf", height=4, width=4)
par(mfrow=c(1,1), mar=c(2,3,2,.5), tcl=-0.25, cex=0.9)
plot(x=yr.st:(n.yrs-1), y=pe[yr.st:(n.yrs-1)], type="n", ylim=c(-1, 1), 
     xaxs="i", yaxs="i", xaxt="n", yaxt="n", xlab="", ylab="")
abline(v=seq(10, 30, by=10), lty=3, col=rgb(0,0,0,.4))
abline(h=c(-0.5, 0, 0.5), lty=3, col=rgb(0,0,0,.4))
abline(v=0, lty=2, col=rgb(0,0,0,.75))
axis(side=1, at=seq(0, 30, by=10), labels=seq(0, 30, by=10), cex.axis=0.8, las=1, mgp=c(0.7, 0.2, 0))
axis(side=2, at=c(-1,-0.5,0,0.5,1), cex.axis=0.8, las=0, mgp=c(1, 0.4, 0))
lines(x=yr.st:(n.yrs-1), y=pe[yr.st:(n.yrs-1)], col=rgb(1,0.55,0,1), lwd=2)
polygon(x=c(yr.st:(n.yrs-1), (n.yrs-1):yr.st), col=rgb(1,0.55,0,.5), border=NA, 
        y=c(u95[yr.st:(n.yrs-1)], l95[(n.yrs-1):yr.st]))
mtext(side=1, line=1.1, "Years", cex=0.9)
mtext(side=2, line=1.6, "Immigration mood", cex=0.9)
text(x=0.5, y=0.8, "Net migration increases by 2 SD in year 0", cex=0.85, pos=4)
mtext(text="Existing immigrant stock is average", side=3, line=0.4, cex=0.9, adj=0)
dev.off()


## Figure 3c

mod = sal.mods[[1]]
dem.mod = imm.stck.mod.1
round(coef(mod), 2)
n.coef =  length(coef(mod))
n.yrs = 30
n.burn = 0

chg.nm = 2 * sd(migrat$Net_migr_within, na.rm=TRUE)
pred.data = data.frame(ImmSal_m1=0, ImmSal_m2=0, 
                       GDPgr=0, Unemp=0, 
                       diff_NM=0,
                       NM=c(rep(0, n.burn), rep(chg.nm, n.yrs-n.burn)), 
                       Noncit_stck_pc=0,
                       FRseat=0, BMIi=0, BMIe=0,
                       DiffFlowsxStock=NA,
                       FlowsxStock=NA,
                       ChgSal=NA, ImmSal=NA
)
n.sims = 1e3
sim.data = array(as.matrix(pred.data), dim=c(dim(pred.data), n.sims))
mod.sims = mvrnorm(n=n.sims, mu=coef(mod), Sigma=plm::vcovSCC(mod, cluster="time"), empirical=FALSE)
dem.mod.sims = mvrnorm(n=n.sims, mu=coef(dem.mod), Sigma=plm::vcovSCC(dem.mod, cluster="time"), empirical=FALSE)
dem.hat = matrix(NA, ncol=n.yrs, nrow=n.sims)
yhat = matrix(NA, ncol=n.yrs, nrow=n.sims)
sup.hat = matrix(NA, ncol=n.yrs, nrow=n.sims)

for(i in 1:(n.yrs-1)) {
  for(k in 1:n.sims) {
    dem.hat[k, i] = dem.mod.sims[k,] %*% as.matrix(sim.data[i, c(7,6), k])
    sim.data[i+1, match("Noncit_stck_pc", names(pred.data)), k] = dem.hat[k, i] + rnorm(1, 0, sigma(dem.mod))
    sim.data[i, match("FlowsxStock", names(pred.data)), k] = sim.data[i, match("NM", names(pred.data)), k] * 
      sim.data[i, match("Noncit_stck_pc", names(pred.data)), k]
    sim.data[i+1, match("diff_NM", names(pred.data)), k] = sim.data[i+1, match("NM", names(pred.data)), k] - 
      sim.data[i, match("NM", names(pred.data)), k]
    sim.data[i, match("DiffFlowsxStock", names(pred.data)), k] = sim.data[i, match("diff_NM", names(pred.data)), k] * 
      sim.data[i, match("Noncit_stck_pc", names(pred.data)), k]
    sim.data[i, match("ChgSal", names(pred.data)), k] = sum(sim.data[i, 1:n.coef, k] * coef(mod))
    sim.data[i, match("ImmSal", names(pred.data)), k] = sim.data[i, match("ImmSal_m1", names(pred.data)), k] + 
      sim.data[i, match("ChgSal", names(pred.data)), k]
    sim.data[i+1, match("ImmSal_m1", names(pred.data)), k] = sim.data[i, match("ImmSal", names(pred.data)), k]
    sim.data[i+1, match("ImmSal_m2", names(pred.data)), k] = sim.data[i, match("ImmSal_m1", names(pred.data)), k]    
    sim.data[i, match("ChgSal", names(pred.data)), k] = mod.sims[k,] %*% (as.matrix(sim.data[i, 1:n.coef, k]))
    sup.hat[k, i] = sim.data[i, match("ImmSal_m1", names(pred.data)), k] + sim.data[i,match("ChgSal", names(pred.data)), k]
  }
}

pe = apply(sup.hat, 2, mean, na.rm=TRUE)
pe.sd = apply(sup.hat, 2, sd, na.rm=TRUE)
u95 = pe + 1.96 * pe.sd
l95 = pe - 1.96 * pe.sd

yr.st = 1
yr1.off = pe[yr.st]
pe = pe-yr1.off
u95 = u95-yr1.off
l95 = l95-yr1.off

pdf("fig3c.pdf", height=4, width=4)
par(mfrow=c(1,1), mar=c(2,3,2,.5), tcl=-0.25, cex=0.9)
plot(x=yr.st:(n.yrs-1), y=pe[yr.st:(n.yrs-1)], type="n", ylim=c(-1, 1), 
     xaxs="i", yaxs="i", xaxt="n", yaxt="n", xlab="", ylab="")
abline(v=seq(10, 30, by=10), lty=3, col=rgb(0,0,0,.4))
abline(h=c(-0.5, 0, 0.5), lty=3, col=rgb(0,0,0,.4))
abline(v=0, lty=2, col=rgb(0,0,0,.75))
axis(side=1, at=seq(0, 30, by=10), labels=seq(0, 30, by=10), cex.axis=0.8, las=1, mgp=c(0.7, 0.2, 0))
axis(side=2, at=c(-1,-0.5,0,0.5,1), cex.axis=0.8, las=0, mgp=c(1, 0.4, 0))
lines(x=yr.st:(n.yrs-1), y=pe[yr.st:(n.yrs-1)], col=rgb(1,0.55,0,1), lwd=2)
polygon(x=c(yr.st:(n.yrs-1), (n.yrs-1):yr.st), col=rgb(1,0.55,0,.5), border=NA, 
        y=c(u95[yr.st:(n.yrs-1)], l95[(n.yrs-1):yr.st]))
mtext(side=1, line=1.1, "Years", cex=0.9)
mtext(side=2, line=1.6, "Immigration concern", cex=0.9)
text(x=0.5, y=0.8, "Net migration increases by 2 SD in year 0", cex=0.85, pos=4)
mtext(text="Existing immigrant stock is average", side=3, line=0.4, cex=0.9, adj=0)
dev.off()


## Figure 2d

mod = perc.mods[[2]]
dem.mod = imm.stck.mod.1
round(coef(mod), 2)
n.coef =  length(coef(mod))
n.yrs = 30
n.burn = 0

chg.immin = 2 * sd(migrat$Immig_pc_within, na.rm=TRUE)
pred.data = data.frame(ImmSup_m1=0, ImmSup_m2=0, 
                       GDPgr=0, Unemp=0, 
                       diff_ImmIn=0,
                       ImmIn=c(rep(0, n.burn), rep(chg.immin, n.yrs-n.burn)), 
                       Noncit_stck_pc=0,
                       FRseat=0, BMIi=0, BMIe=0,
                       DiffFlowsxStock=NA,
                       FlowsxStock=NA,
                       ChgSup=NA, ImmSup=NA
)
n.sims = 1e3
sim.data = array(as.matrix(pred.data), dim=c(dim(pred.data), n.sims))
mod.sims = mvrnorm(n=n.sims, mu=coef(mod), Sigma=plm::vcovSCC(mod, cluster="time"), empirical=FALSE)
dem.mod.sims = mvrnorm(n=n.sims, mu=coef(dem.mod), Sigma=plm::vcovSCC(dem.mod, cluster="time"), empirical=FALSE)
dem.hat = matrix(NA, ncol=n.yrs, nrow=n.sims)
yhat = matrix(NA, ncol=n.yrs, nrow=n.sims)
sup.hat = matrix(NA, ncol=n.yrs, nrow=n.sims)

for(i in 1:(n.yrs-1)) {
  for(k in 1:n.sims) {
    dem.hat[k, i] = dem.mod.sims[k,] %*% as.matrix(sim.data[i, c(7,6), k])
    sim.data[i+1, match("Noncit_stck_pc", names(pred.data)), k] = dem.hat[k, i] + rnorm(1, 0, sigma(dem.mod))
    sim.data[i, match("FlowsxStock", names(pred.data)), k] = sim.data[i, match("ImmIn", names(pred.data)), k] * 
      sim.data[i, match("Noncit_stck_pc", names(pred.data)), k]
    sim.data[i+1, match("diff_ImmIn", names(pred.data)), k] = sim.data[i+1, match("ImmIn", names(pred.data)), k] - 
      sim.data[i, match("ImmIn", names(pred.data)), k]
    sim.data[i, match("DiffFlowsxStock", names(pred.data)), k] = sim.data[i, match("diff_ImmIn", names(pred.data)), k] * 
      sim.data[i, match("Noncit_stck_pc", names(pred.data)), k]
    sim.data[i, match("ChgSup", names(pred.data)), k] = sum(sim.data[i, 1:n.coef, k] * coef(mod))
    sim.data[i, match("ImmSup", names(pred.data)), k] = sim.data[i, match("ImmSup_m1", names(pred.data)), k] + 
      sim.data[i, match("ChgSup", names(pred.data)), k]
    sim.data[i+1, match("ImmSup_m1", names(pred.data)), k] = sim.data[i, match("ImmSup", names(pred.data)), k]
    sim.data[i+1, match("ImmSup_m2", names(pred.data)), k] = sim.data[i, match("ImmSup_m1", names(pred.data)), k]    
    sim.data[i, match("ChgSup", names(pred.data)), k] = mod.sims[k,] %*% (as.matrix(sim.data[i, 1:n.coef, k]))
    sup.hat[k, i] = sim.data[i, match("ImmSup_m1", names(pred.data)), k] + sim.data[i,match("ChgSup", names(pred.data)), k]
  }
}

pe = apply(sup.hat, 2, mean, na.rm=TRUE)
pe.sd = apply(sup.hat, 2, sd, na.rm=TRUE)
u95 = pe + 1.96 * pe.sd
l95 = pe - 1.96 * pe.sd
yr.st = 1
yr1.off = pe[yr.st]
pe = pe-yr1.off
u95 = u95-yr1.off
l95 = l95-yr1.off

pdf("fig2d.pdf", height=4, width=4)
par(mfrow=c(1,1), mar=c(2,3,2,.5), tcl=-0.25, cex=0.9)
plot(x=yr.st:(n.yrs-1), y=pe[yr.st:(n.yrs-1)], type="n", ylim=c(-1, 1), 
     xaxs="i", yaxs="i", xaxt="n", yaxt="n", xlab="", ylab="")
abline(v=seq(10, 30, by=10), lty=3, col=rgb(0,0,0,.4))
abline(h=c(-0.5, 0, 0.5), lty=3, col=rgb(0,0,0,.4))
abline(v=0, lty=2, col=rgb(0,0,0,.75))
axis(side=1, at=seq(0, 30, by=10), labels=seq(0, 30, by=10), cex.axis=0.8, las=1, mgp=c(0.7, 0.2, 0))
axis(side=2, at=c(-1,-0.5,0,0.5,1), cex.axis=0.8, las=0, mgp=c(1, 0.4, 0))
lines(x=yr.st:(n.yrs-1), y=pe[yr.st:(n.yrs-1)], col=rgb(1,0.55,0,1), lwd=2)
polygon(x=c(yr.st:(n.yrs-1), (n.yrs-1):yr.st), col=rgb(1,0.55,0,.5), border=NA, 
        y=c(u95[yr.st:(n.yrs-1)], l95[(n.yrs-1):yr.st]))
mtext(side=1, line=1.1, "Years", cex=0.9)
mtext(side=2, line=1.6, "Immigration mood", cex=0.9)
text(x=0.5, y=0.8, "Immigration inflows increase by 2 SD in year 0", cex=0.85, pos=4)
mtext(text="Existing immigrant stock is average", side=3, line=0.4, cex=0.9, adj=0)
dev.off()


## Figure 3d

mod = sal.mods[[2]]
dem.mod = imm.stck.mod.1
round(coef(mod), 2)
n.coef =  length(coef(mod))
n.yrs = 30
n.burn = 0

chg.immin = 2 * sd(migrat$Immig_pc_within, na.rm=TRUE)
pred.data = data.frame(ImmSal_m1=0, ImmSal_m2=0, 
                       GDPgr=0, Unemp=0, 
                       diff_ImmIn=0,
                       ImmIn=c(rep(0, n.burn), rep(chg.immin, n.yrs-n.burn)), 
                       Noncit_stck_pc=0,
                       FRseat=0, BMIi=0, BMIe=0,
                       DiffFlowsxStock=NA,
                       FlowsxStock=NA,
                       ChgSal=NA, ImmSal=NA
)
n.sims = 1e3
sim.data = array(as.matrix(pred.data), dim=c(dim(pred.data), n.sims))
mod.sims = mvrnorm(n=n.sims, mu=coef(mod), Sigma=plm::vcovSCC(mod, cluster="time"), empirical=FALSE)
dem.mod.sims = mvrnorm(n=n.sims, mu=coef(dem.mod), Sigma=plm::vcovSCC(dem.mod, cluster="time"), empirical=FALSE)
dem.hat = matrix(NA, ncol=n.yrs, nrow=n.sims)
yhat = matrix(NA, ncol=n.yrs, nrow=n.sims)
sup.hat = matrix(NA, ncol=n.yrs, nrow=n.sims)

for(i in 1:(n.yrs-1)) {
  for(k in 1:n.sims) {
    dem.hat[k, i] = dem.mod.sims[k,] %*% as.matrix(sim.data[i, c(7,6), k])
    sim.data[i+1, match("Noncit_stck_pc", names(pred.data)), k] = dem.hat[k, i] + rnorm(1, 0, sigma(dem.mod))
    sim.data[i, match("FlowsxStock", names(pred.data)), k] = sim.data[i, match("ImmIn", names(pred.data)), k] * 
      sim.data[i, match("Noncit_stck_pc", names(pred.data)), k]
    sim.data[i+1, match("diff_ImmIn", names(pred.data)), k] = sim.data[i+1, match("ImmIn", names(pred.data)), k] - 
      sim.data[i, match("ImmIn", names(pred.data)), k]
    sim.data[i, match("DiffFlowsxStock", names(pred.data)), k] = sim.data[i, match("diff_ImmIn", names(pred.data)), k] * 
      sim.data[i, match("Noncit_stck_pc", names(pred.data)), k]
    sim.data[i, match("ChgSal", names(pred.data)), k] = sum(sim.data[i, 1:n.coef, k] * coef(mod))
    sim.data[i, match("ImmSal", names(pred.data)), k] = sim.data[i, match("ImmSal_m1", names(pred.data)), k] + 
      sim.data[i, match("ChgSal", names(pred.data)), k]
    sim.data[i+1, match("ImmSal_m1", names(pred.data)), k] = sim.data[i, match("ImmSal", names(pred.data)), k]
    sim.data[i+1, match("ImmSal_m2", names(pred.data)), k] = sim.data[i, match("ImmSal_m1", names(pred.data)), k]    
    sim.data[i, match("ChgSal", names(pred.data)), k] = mod.sims[k,] %*% (as.matrix(sim.data[i, 1:n.coef, k]))
    sup.hat[k, i] = sim.data[i, match("ImmSal_m1", names(pred.data)), k] + sim.data[i,match("ChgSal", names(pred.data)), k]
  }
}

pe = apply(sup.hat, 2, mean, na.rm=TRUE)
pe.sd = apply(sup.hat, 2, sd, na.rm=TRUE)
u95 = pe + 1.96 * pe.sd
l95 = pe - 1.96 * pe.sd
yr.st = 1
yr1.off = pe[yr.st]
pe = pe-yr1.off
u95 = u95-yr1.off
l95 = l95-yr1.off

pdf("fig3d.pdf", height=4, width=4)
par(mfrow=c(1,1), mar=c(2,3,2,.5), tcl=-0.25, cex=0.9)
plot(x=yr.st:(n.yrs-1), y=pe[yr.st:(n.yrs-1)], type="n", ylim=c(-1, 1), 
     xaxs="i", yaxs="i", xaxt="n", yaxt="n", xlab="", ylab="")
abline(v=seq(10, 30, by=10), lty=3, col=rgb(0,0,0,.4))
abline(h=c(-0.5, 0, 0.5), lty=3, col=rgb(0,0,0,.4))
abline(v=0, lty=2, col=rgb(0,0,0,.75))
axis(side=1, at=seq(0, 30, by=10), labels=seq(0, 30, by=10), cex.axis=0.8, las=1, mgp=c(0.7, 0.2, 0))
axis(side=2, at=c(-1,-0.5,0,0.5,1), cex.axis=0.8, las=0, mgp=c(1, 0.4, 0))
lines(x=yr.st:(n.yrs-1), y=pe[yr.st:(n.yrs-1)], col=rgb(1,0.55,0,1), lwd=2)
polygon(x=c(yr.st:(n.yrs-1), (n.yrs-1):yr.st), col=rgb(1,0.55,0,0.5), border=NA, 
        y=c(u95[yr.st:(n.yrs-1)], l95[(n.yrs-1):yr.st]))
mtext(side=1, line=1.1, "Years", cex=0.9)
mtext(side=2, line=1.6, "Immigration concern", cex=0.9)
text(x=0.5, y=0.8, "Immigration inflows increase by 2 SD in year 0", cex=0.85, pos=4)
mtext(text="Existing immigrant stock is average", side=3, line=0.4, cex=0.9, adj=0)
dev.off()


